/*
 * Voronoicavities.cpp
 *
 *  Created on: May 26, 2012
 *      Author: marchi
 */

#include "Voronoi_cavities.h"

namespace Cavities {

Voronoi_cavities::Voronoi_cavities(int n, const vector<string> & s, bool bH): Voronoi(),hstg(vector<Histogram>(1000)) {
	nr=n;
	Vol=new double [nr];
	Neighs=new vector<int> [nr];
	Surface=new vector<double> [nr];
	label=new string [nr];
	types=new int [nr];

	for(int i=0;i<nr;i++){
		label[i]=s[i];
		if(bH || label[i].compare(0,1,h,0,1)) cindex.push_back(i);
	}
}

void Voronoi_cavities::getDelaunay(isTypes & isWhat){
	int ia=0;
	if(!Mycon) return;
	c_loop_order_periodic vl(*Mycon,*porder);
	voronoicell_neighbor c;
	vector<double> vertx;
	cv.clear();
	cvnei.clear();
	cvnd.clear();
	set<Vertex> cvE;
	if(vl.start())
		do {
			vector<int> nodes;
			if(Mycon->compute_cell(c,vl)) {
				int p=c.p;

				c.vertices(vertx);
				double x,y,z,r;
				int iD;
				vl.pos(iD,x,y,z,r);
				for(int m=0;m<p;m++){
					int xx=m*DIM+XX;
					int yy=xx+1;
					int zz=yy+1;
					double myr=sqrt(vertx[xx]*vertx[xx]+vertx[yy]*vertx[yy]+vertx[zz]*vertx[zz]);
					nodes.push_back(cindex[ia]);
					if(c.nu[m]>0) {
						int j=0;
						while(j<c.nu[m]-1) nodes.push_back(cindex[c.ne[m][j++]]);
						nodes.push_back(cindex[c.ne[m][j]]);
					};

					try{if(nodes.size() != 4) throw string(" Warning Number of nodes is not 4 but ");}
					catch(const string & s){cout << s << nodes.size() << endl;}

					Vertex tm(nodes);
					Vertex tm2(tm,Top.gbG());
					if(myr-r < 0.0) {
						cvE.insert(tm);
						map<const Vertex,CavityData >::iterator it=cv.find(tm);
						if(it != cv.end()) cv.erase(it);
					} else if(cvE.find(tm) == cvE.end() && isWhat(tm2)){
						cv[tm](myr-r,Dvect(vertx[xx]+x,vertx[yy]+y,vertx[zz]+z),tm2);
					}
					nodes.clear();
				}
			}
			ia++;
		} while(vl.inc());

	map<const Vertex,CavityData >::iterator it=cv.begin();
	for(;it != cv.end();++it){
		SetCavityData tmp(*it);
		set<SetCavityData>::iterator its=cvnd.find(tmp);
		if(its != cvnd.end()) {
			SetCavityData tmp0=*its;
			tmp0.Merge(*it,Top.gbG());
			cvnd.erase(its);
			cvnd.insert(tmp0);
		}else cvnd.insert(tmp);
	}
	cout << "----- Identified " << cv.size()- cvnd.size()<< " quasi-degenerate vertices over "<< cv.size() << " -----"<< endl;
	set<SetCavityData>::iterator itc=cvnd.begin();
	for(;itc != cvnd.end();++itc){
		SetCavityData tmp(*itc);
		double R=tmp.gMinRadius();
		vector<int> tm=tmp.gtyId();
		for(unsigned int n=0;n<tm.size();n++) hstg[tm[n]](R);
	}
}
void Voronoi_cavities::WriteIt(std::ofstream & fout){
	vector<vector<int> > tmp=Top.gvG();
	for(int n=0;n< Top.gbGMax();n++){
		double sum=hstg[n].TotHist();
		hstg[n].Avg();
		string ss=" ";
		ss+=Top.AtomResidue(tmp[n][0])+ " : ";
		if(n)
			for(unsigned int m=0;m<tmp[n].size();m++){
				ss+=label[tmp[n][m]]+ " ";
			}
		fout <<"& " << n << "  sum = " << sum << ss << endl;
		fout << hstg[n];
	}
}

void Voronoi_cavities::copyhgtotot(){
	vector<Histogram>::iterator ith=hstg.begin();
	for(;ith != hstg.end();++ith){
		copy(ith->gHG().begin(),ith->gHG().end(),back_inserter(tot_hg));
	}
}

void Voronoi_cavities::copytottohg(){
	vector<Histogram>::iterator ith=hstg.begin();
	vector<double>::iterator it0=tot_hg.begin();
	for(;ith != hstg.end();++ith){
		ith->sHG(it0);
	}
}

Voronoi_cavities::~Voronoi_cavities() {
	// TODO Auto-generated destructor stub
}

} /* namespace Cavities */
